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ABSTRACT 

The  objectives  of  this  project  are 

•  To  investigate  the  limits  of  empirical  matched  field  processing  and  other  coherent  array  detection  and 
parameter  estimation  methods  as  receiver  aperture  size  increases  from  a  few  kilometers  to  many  hundreds 
of  kilometers. 

•  To  investigate  techniques  for  extending  the  geographical  source-region  footprint  over  which  empirical 
matched  field  processing  and  other  coherent  calibrated  methods  apply. 

Under  a  previous  contract,  we  demonstrated  the  ability  of  empirical  matched  field  processing  to  classify  mining 
explosions  by  originating  mine,  using  data  from  a  single  small-aperture  array  (ARCES)  applied  to  mining  events  on 
the  Kola  Peninsula.  In  the  current  contract,  we  have  chosen  central  Asia  as  our  study  region  to  assure  programmatic 
relevance  and  to  exploit  the  large  belts  of  natural  and  man-made  seismicity  required  for  a  test  of  our  processing 
strategy.  Data  from  suitable  networks  available  for  the  study  are  from  the  four  Kazakhstan  arrays  (MKAR,  KKAR, 
ABKAR,  BRVK)  and  the  Kyrgyzstan  network.  Our  overall  approach  will  be  to  start  with  small  apertures  (the 
individual  Kazakhstan  arrays  considered  separately)  to  check  the  reproducibility  of  results  we  obtained  for  the 
European  Arctic,  then  to  extend  the  processing  strategy  first  to  a  medium  aperture  (the  Kyrgyzstan  network;  200  km 
aperture)  and  then  a  large  aperture  (the  four  Kazakhstan  arrays  considered  as  a  single  coherent  aperture;  >1000  km). 
Simultaneously  we  will  investigate  methods  to  expand  the  source  region  footprint  over  which  calibrations  for 
coherent  methods  apply. 

We  have  expanded  our  work  on  the  dataset  of  mining  explosions  in  the  Kola  Peninsula  to  investigate  further  the 
potential  of  the  method  and  to  validate  the  preliminary  results  obtained  under  a  previous  contract  (BAA03-27).  This 
validation  has  consisted  of  analyzing  several  hundred  seismic  waveforms  from  available  local  stations  (Apatity  and 
Lovozero)  to  make  additional  checks  on  the  original  ground  truth  information  provided  by  the  mining  authorities. 
Cross  validation  also  was  introduced  to  avoid  circularity  in  the  data  analysis.  Through  the  use  of  empirical  matched 
field  processing  we  were  able  correctly  to  classify  539  of  549  explosions  by  originating  mine  among  10  closely 
spaced  mines,  using  short  segments  (<10  seconds,  filtered  in  the  2.5-12.5  Hz  band)  of  Pn  waves  observed  at  the 
ARCES  array.  The  mines  are  separated  by  as  little  as  3  kilometers  and  observed  at  a  340-410  kilometer  range.  Such 
classification  performance  is  far  superior  to  that  obtained  by  traditional  methods.  Matched  field  processing  should 
find  application  in  screening  mining  explosions  in  Comprehensive  Nuclear-Test-Ban  Treaty  monitoring  applications 
and  in  studies  of  natural  seismicity. 
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OBJECTIVES 

The  objectives  of  this  project  are 

•  To  investigate  the  limits  of  empirical  matched  field  processing  and  other  coherent  array  detection  and 
parameter  estimation  methods  as  receiver  aperture  size  increases  from  a  few  kilometers  to  many  hundreds 
of  kilometers. 

•  To  investigate  techniques  for  extending  the  geographical  source-region  footprint  over  which  empirical 
matched  field  processing  and  other  coherent  calibrated  methods  apply. 

We  will  begin  by  reanalyzing  data  from  the  European  Arctic  in  order  to  reconfirm  the  potential  of  empirical 
matched  field  processing  that  has  been  indicated  by  our  work  under  a  previous  contract.  We  then  will  proceed  to 
study  the  central  Asia  region  to  assure  programmatic  relevance  and  to  exploit  the  large  belts  of  natural  and 
man-made  seismicity  to  test  the  applicability  of  the  technique  to  diffuse  seismicity. 

RESEARCH  ACCOMPLISHED 

Introduction 

Current  operational  seismic  array  processing  methods  for  detection  (beamforming)  and  parameter  estimation 
(frequency- wavenumber  analysis)  have  changed  little  since  their  introduction  in  the  1960s  and  1970s.  These 
methods  rely  upon  a  plane-wave  assumption  for  predicting  the  spatial  amplitude  and  phase  structure  of  seismic 
waves  incident  upon  an  array  aperture.  Scattering  and  refraction  in  strongly  heterogeneous  seismic  propagation 
media  constrain  the  size  of  usable  coherent  processing  apertures  under  the  plane-wave  assumption  to  a  few 
wavelengths.  This  constraint  severely  limits  the  spatial  resolution  of  beamforming  and  frequency-wavenumber 
(F-K)  methods.  In  addition,  the  spatial  correlation  of  ambient  seismic  noise  constrains  the  minimum  usable  element 
separation.  The  combination  of  these  constraints  bounds  the  maximum  number  of  usable  sensors  and  places  a 
fundamental  limit  on  coherent  processing  gain  through  beamforming. 

Existing  beamforming  and  F-K  estimation  methods  make  some  attempt  to  compensate  for  the  non-ideal  spatial 
structure  of  seismic  waves.  For  example,  it  is  common  to  apply  an  amplitude  correction  to  estimated  spatial 
covariance  matrices  when  using  the  indirect  approach  to  F-K  spectrum  estimation  (i.e.,  by  estimating  the  the  spatial 
covariance  function  first,  then  evaluating  its  Fourier  transform  to  obtain  a  wavenumber  spectrum).  Each  element  of 
the  matrix  is  normalized  by  the  square  root  of  the  diagonal  elements  in  the  same  row  and  column.  In  direct  methods 
(e.g.,  Kvaema  and  Ringdal,  1986),  it  is  now  standard  to  integrate  the  F-K  spectrum  over  a  band  of  frequencies  to 
stabilize  the  narrowband  estimates  that  individually  have  high  variance.  That  variance  is  controlled,  in  part,  by  the 
very  small  time-bandwidth  product  of  signals  in  short  analysis  windows.  However,  it  is  likely  that  variance  in  the 
F-K  estimate  from  frequency  to  frequency  also  is  a  function  of  wavefield  scattering. 

Recent  advances  in  standardization  of  analysis  windows  and  frequency  bands  for  particular  source  regions  have 
produced  dramatic  reductions  in  the  variance  of  azimuth/slowness  estimates  derived  from  F-K  spectra  (Gibbons  et 
al.,  2005),  at  least  in  frequency  bands  (below  6-8  Hz)  where  the  wavefield  is  approximately  coherent  across 
regional-array  apertures.  Finally,  it  is  common  to  apply  post-processing  (i.e.,  post-F-K)  vector  slowness  corrections 
(Schweitzer,  2001)  to  azimuth/slowness  measurements  in  an  attempt  to  remove  frequency-dependent  biases  caused 
(presumably)  by  unmodeled  wavefield  refraction. 

Recently,  detection  algorithms  that  exploit  aperture-level  calibrations  have  been  developed  that  often  substantially 
reduce  detection  thresholds  in  beamforming  operations.  Correlation  detectors  (Gibbons  and  Ringdal,  2006)  and 
subspace  detectors  (Harris,  2006)  use  previously  observed  waveforms  from  events  at  specific  sources  in  sensitive 
algorithms  to  detect  subsequent  events  at  those  same  sources.  These  algorithms  rely  upon  the  fact  that  the  temporal 
and  spatial  structures  of  signals  from  these  sources  often  are  repeatable;  past  events  constitute  a  spatio-temporal 
calibration  across  an  array  or  network  of  sensors  for  future  events.  Such  detectors  simultaneously  detect,  locate,  and 
identify  events  as  being  consistent  with  previous  events  at  the  same  source  location,  with  the  same  source  time 
history  and  mechanism. 
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However,  while  the  full  exploitation  of  the  spatial  and  temporal  structure  of  the  signal  frequently  leads  to  detectors 
one  to  two  orders  of  magnitude  more  sensitive  than  simple  power  detectors  (with  conventional  beamforming),  such 
exploitation  also  presents  two  barriers  to  widespread  application.  Correlation  methods  generally  have  a 
source-region  geographic  footprint  one  to  two  wavelengths  across  (Harris,  1991),  which  can  be  ameliorated  to  some 
extent  with  the  more  general  waveform  representation  of  a  subspace  detector  (Harris,  2006).  Correlation  methods 
correct  incoherence  across  the  receiver  aperture,  but  trade  it  for  incoherence  across  the  source  region;  the  reciprocity 
principle  is  at  play  in  the  correlation  detection  strategy.  Correlation  methods  also  are  sensitive  to  source  mechanism 
and  time  history,  and  can  fail  when  these  attributes  differ  from  event  to  event. 

The  sensitivities  of  F-K  methods  and  conventional  beamforming  to  the  plane  wave  assumption  and  of  correlation 
methods  to  source  variability  leads  us  to  consider  another  processing  strategy  based  on  temporally-incoherent,  but 
spatially  coherent  signal  processing  extended  with  subspace  techniques.  Our  general  approach  is  to  adapt 
narrowband  empirical  matched  field  processing,  originally  conceived  for  underwater  acoustic  applications 
(Baggeroer  et  al.,  1993;  Fialkowski  et  al.,  2000)  to  the  seismic  monitoring  task.  Our  processing  agenda  consists  of 
( 1 )  breaking  the  observed  signal  into  narrow  bands,  (2)  using  empirical  matched  field  methods  to  combine  the 
near-monochromatic  resultant  waveforms  coherently  across  an  aperture  to  achieve  processing  gain,  and 
(3)  integrating  the  resulting  power  incoherently  over  the  narrow  bands  to  achieve  a  wideband  result.  We  anticipate 
that  this  approach  will  ameliorate  variations  in  source  time  history  from  event  to  event  that  defeat  waveform 
correlation  methods. 

Matched  field  processing  derives  its  name  from  the  realization  that  conventional  narrowband  beamforming  is  a 
spatial  matched  filtering  operation  for  plane  waves,  and  that  more  general  forms  of  spatial  matched  filtering  are 
possible.  In  underwater  acoustics,  the  velocity  structure  of  the  medium  is  far  less  variable  than  the  seismic 
propagation  medium  and  frequently  known  in  detail.  In  such  circumstances,  it  is  possible  to  compute  the  detailed 
phase  and  amplitude  structure  of  a  narrowband  wavefield  incident  on  an  array  and  apply  the  computed  structure  as 
the  focusing  kernel  in  a  beamforming  operation.  This  approach  works  well  in  the  SOund  Fixing  And  Ranging 
(SOFAR)  channel  waveguide  and  allows  coherent  processing  across  arrays  to  be  expanded  beyond  the  physical 
aperture  limits  implied  by  the  assumption  of  a  planar  structure  (locally  the  field  is  planar).  A  clever  empirical 
approach  to  estimating  the  structure  of  the  narrowband  acoustic  wavefield  directly  from  observed  wavefields  was 
attempted  (Fialkowski  et  al.,  2000)  but  largely  failed  because  of  the  spatial  nonstationarity  of  the  source. 

Reanalyzing  Mining  Data  from  the  European  Arctic 

Under  a  previous  contract,  DE-FC52-03NA99517  (Kvaema  et  al.,  2007),  we  demonstrated  the  ability  of  empirical 
matched  field  processing  to  classify  mining  explosions  by  originating  mine,  using  data  from  a  single  small-aperture 
array  (ARCES)  applied  to  mining  events  on  the  Kola  Peninsula. 

We  have  continued  our  work  on  this  data  set  for  the  mining  explosions  at  the  Kola  Peninsula  in  order  to  further 
validate  the  preliminary  results  summarized  by  Kvaema  et  al.  (2007).  This  validation  has  consisted  of  analyzing  over 
500  seismic  waveforms  from  available  local  stations  (Apatity  and  Lovozero)  to  make  additional  checks  on  the 
original  ground  truth  information  provided  by  the  mining  authorities.  We  have  also  carried  out  classification  tests 
using  cross-validation  to  avoid  issues  of  circularity. 

We  study  the  European  arctic  region  (Figure  1)  where  the  ARCES  array,  located  in  northern  Norway,  is  340-410 
kilometers  from  two  groups  of  mines  in  the  Kola  peninsula,  Russia:  the  Olenegorsk  (01-05)  and  Khibiny  (K1-K5) 
groups.  In  monitoring  for  signals  from  distant  nuclear  events,  the  array  observes  thousands  of  mining  explosions 
annually,  which  can  be  screened  if  attributed  with  a  high  level  of  confidence  to  their  originating  mines.  However, 
the  array,  with  a  3-kilometer  aperture,  is  too  small  to  resolve  the  individual  mines.  The  array  resolution  (Rayleigh 
criterion)  is  determined  by  the  separation  of  half-power  (3  dB)  points  of  the  main  lobe  of  the  array  wavenumber 
response.  In  Figure  I ,  the  3  dB  points  of  the  array  response  (at  4,  8,  and  12  Hz)  for  first-arriving  Pn  waves  are 
projected  onto  the  array’s  geographic  field  of  view;  the  array  has  been  steered  to  one  particular  mine  (K2).  All  10 
mines  fall  within  the  3  dB  contours  even  at  12  Hz. 
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Figure  1.  The  mines  of  the  Khibiny  and  Olenegorsk  regions  are  too  close  for  explosions  to  be  attributed  to 

specific  mines  on  the  basis  of  wavenumber  (F-K)  spectrum  direction  estimates  made  from  ARCES 
array  (upper  left)  observations  of  Pn  waves  (Figure  2).  ARCES  resolution  is  indicated  on  the  map 
at  lower  left,  which  shows  the  half-power  contours  of  the  array  response  at  three  frequencies  when 
the  array  is  steered  to  the  Rasvumchorr  mine  (K2). 

We  consider  separating  explosions  from  the  10  mines  using  only  the  first-arriving  Pn  wave.  Figure  2  displays 
ARCES  center  station  observations  of  37  explosions  originating  at  the  Norpakh  mine  (K5).  Note  the  large  degree  of 
signal  variation,  which  complicates  the  use  of  correlation  methods.  For  our  analysis,  we  selected  549  explosions 
attributed  to  specific  mines  using  reports  from  the  mine  operators,  validated  with  data  from  stations  (APA,  LVZ) 
local  to  the  mines. 


The  classical  technique  for  estimating  incident  wave  velocity  and  direction  is  the  F-K  spectrum, 

2 


P{o),k)  = 


=  R{o))e 


(1) 


which  maps  the  energy  incident  upon  the  array  as  a  function  of  frequency  CO  and  wavenumber  (direction)  k  by 
computing  a  multidimensional  Fourier  transform  of  the  signals  recorded  at  the  sensor  locations  x  ^ .  The 

F-K  spectrum  can  be  written  as  the  indicated  quadratic  form  between  the  steering  vector 
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which  expresses  the  spatial  phase  and  amplitude  structure  of  incident 


monochromatic  plane  waves,  and  the  spatial  covariance  matrix  R  ,  containing  the  cross-spectral  information 
^mn  (^)  “  ( dt )( )e'^ dt)  among  the  sensor  waveforms. 


Figure  2.  The  ensemble  of  events  used  to  calibrate  the  Norpakh  mine  (K5).  Note  that  the 

first-arriving  P  waves  (see  inset)  show  a  degree  of  variation  that  makes  classification  by 
waveform  correlation  methods  infeasible. 


It  is  common  to  integrate  the  energy  over  frequency  to  produce  a  map  of  total  incident  energy  as  a  function  of  vector 
slowness,  s  , 
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P(?)=  ^P{co,cdsyico 

This  function  aggregates  energy  incoherently  across  frequency,  which  is  the  basis  for  wideband  estimation 
algorithms  not  sensitive  to  the  specific  temporal  structure  of  signals. 

The  direction  to  the  source  of  an  observed  signal  is  estimated  by  maximizing  ^(?)  over  vector  slowness.  In 
practice,  the  wavefield  often  is  strongly  refracted  leading  to  significant  bias  in  direction  estimates.  It  is  common  to 
compensate  refraction  effects  with  empirical  corrections  to  the  estimated  slowness  derived  from  events  of  known 
location  (Schweitzer,  2001).  This  calibration  approach  effectively  replaces  the  theoretical  steering  vector  for  a 
particular  direction  and  velocity  with  a  corrected  vector,  although  one  that  still  conforms  to  a  plane  wave  signal. 

Improvements  in  estimation  accuracy  result  from  improving  the  representation  of  the  incident  wavefield  structure, 
embodied  in  the  steering  vectors.  In  underwater  sound  localization  problems,  matched  field  processing  (Baggeroer 
et  al.,  1993)  maximizes  the  quadratic  form  of  Equation  (1)  over  a  manifold  of  steering  vectors  computed  from 
detailed  models  of  the  oceanic  sound  velocity  profile  (SVP).  Sufficiently  accurate  SVP  models  are  not  always 
available,  which  motivated  attempts  to  develop  steering  vectors  empirically  (Fialkowski  et  al.,  2000).  In 
hydroacoustic  applications,  signal  sources  frequently  are  ships  producing  continuous  waveforms  modeled  as 
narrowband  random  processes.  When  a  single  source  is  present,  the  steering  vectors  are  obtained  as  the  principal 
eigenvectors  of  sample  covariance  matrices  estimated  from  array  observations  of  the  continuous  wavefield. 

In  adapting  the  matched  field  processing  approach  to  seismic  applications,  an  empirical  approach  is  required,  since 
geophysical  models  of  the  propagation  environment  (the  Earth)  are  not  sufficiently  detailed  or  accurate  to  predict 
wavefield  structure,  except  at  very  low  frequencies  (<  ~0.1  Hz).  We  are  interested  in  improving  performance  at 
frequencies  as  high  as  12  Hz.  Since  the  typical  seismogram  consists  of  many  different  waves  propagating  in  a 
multipath  environment,  we  limit  the  observation  window  to  a  small  portion  of  the  waveform  (approximately  the  first 
ten  seconds)  where  we  can  be  reasonably  sure  that  a  single  wave  is  present.  The  short  window  raises  the  question  of 
how  to  obtain  suitably  stable  estimates  of  spectral  covariance  matrices.  We  solve  this  problem  by  assembling 
ensembles  of  events,  such  as  that  shown  in  Figure  2,  from  repeating  sources  (mines).  The  covariance  matrices  are 
computed  as  the  ensemble  average. 

To  provide  a  baseline  for  array  performance  consistent  with  the  classical  resolution  limit,  we  classify  the  549 
explosions  by  originating  mine  first  by  maximizing  P(?)  over  the  10  plane-wave  steering  vectors  appropriate  to  the 
great-circle  azimuths  to  the  mines.  We  refer  to  this  approach  as  theoretical  plane-wave  (F-K)  classification.  We  also 
perform  the  classification  with  calibrated  plane-wave  steering  vectors,  for  which  the  directions  and  velocities  have 
been  chosen  to  maximize  the  ensemble  power  for  each  mine.  We  refer  to  this  approach  as  empirical  plane  wave 
classification.  Finally,  we  perform  matched  field  classification  using  steering  vectors  obtained  as  the  principal 
eigenvectors  of  the  ensemble  covariance  matrices.  In  all  cases,  we  estimate  the  energy  in  the  2.5-12.5  Hz 
frequency  band,  discretizing  this  band  into  33  subbands  of  0.3125-Hz  width  for  the  covariance  calculations.  We  use 
a  cross-validation  approach  to  avoid  circularity  in  the  estimation  process. 

The  three  algorithms  represent  different  choices  for  representing  the  narrowband  spatial  structure  of  the  incident 
wavefield  under  the  10  source  hypotheses.  To  assess  the  quality  of  the  representations,  we  measure  the  normalized 
energy 
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which  ranges  between  0  and  1  and  is  1  only  if  the  steering  vectors  represent  100%  of  the  signal  energy  incident  on 
the  array.  It  is  the  fraction  of  energy  captured  by  the  specific  representation. 


460 


2008  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


Energy  capture  for  one  mine,  Olenegorsk,  under  the  10  source  hypotheses  is  shown  in  Figure  3.  Both  plane-wave 
classification  methods  have  energy  capture  distributions  centered  about  0.2-0.3  for  Olenegorsk  events  under  the 
Olenegorsk  origin  hypothesis.  The  energy  capture  for  the  alternative  hypotheses,  which  ideally  should  be  0,  is  lower 
for  the  empirical  plane  wave  classifier  than  its  theoretical  counterpart,  suggesting  better  classification  performance 
for  the  calibrated  classifier.  Distributions  for  the  matched  field  method  show  a  very  wide  spread  between  the 
Olenegorsk  hypothesis,  for  which  energy  capture  averages  about  0.8,  and  alternate  hypotheses,  predicting  a  high 
classification  success  rate. 


Figure  3.  Distributions  of  the  matched  field  processing  classification  statistic  (right)  for  the  Olenegorsk  (02) 
mine  population  of  52  events  under  the  10  hypotheses  about  originating  source  (indicated  at  left) 
separate  02  unambiguously  from  the  other  mines.  Frequency  distributions  for  the  theoretical 
plane-wave  F-K  classifier  (left)  are  ambiguous.  Distributions  for  the  empirical  plane  wave  F-K 
classifier  show  slightly  improved  separation  between  the  two  mining  groups. 

Histograms  of  the  relative  frequency  of  classification  for  events  from  each  mine  obtained  by  maximizing  P(?) 
under  the  10  hypotheses  are  shown  in  Figure  4.  In  the  theoretical  plane-wave  case,  a  single  mine,  Rasvumchorr 
(K2),  is  assigned  most  of  the  events.  Observed  Pn  azimuths  generally  are  biased  clockwise  (to  the  south)  for  most 
events  in  this  region.  The  Rasvumchorr  mine  is  the  southernmost  mine,  consequently  the  theoretical  steering  vector 
for  this  mine  presents  the  best  fit  to  the  refracted  wavefields.  The  empirical  plane-wave  classifier  performs 
considerably  better,  generally  assigning  Olenegorsk  group  events  to  Olenegorsk  mines  and  KJiibiny  events  to 
Khibiny  mines,  probably  as  a  consequence  of  calibrating  wideband  refraction  effects.  However,  the  empirical 
plane-wave  classifier  has  large  error  rates  when  distinguishing  among  the  5  mines  in  each  group.  By  contrast  the 
matched-field  classifier  achieves  nearly  complete  separation  of  the  mines  ( 1 .8%  misclassification  rate). 
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